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ABSTRACT 



A bosonic dark matter model is examined in detail via high-resolution simula- 
tions. These bosons have particle mass of order 10~ 22 eV and are non-interacting. 
If they do exist and can account for structure formation, these bosons must be 
condensed into the Bose-Einstein state and described by a coherent wave func- 
tion. This matter, also known as Fuzzy Dark Matter (IHu. Barkana fc Gruzinov 



2000l ). is speculated to be able, first, to eliminate the sub-galactic halos to solve 



the problem of over- abundance of dwarf galaxies, and, second, to produce flat 
halo cores in galaxies suggested by some observations. We investigate this model 
with simulations up to 1024 3 resolution in an 1 h^Mpc box that maintains the 
background matter density Q m = 0.3 and = 0.7. Our results show that the 
extremely light bosonic dark matter (ELBDM) can indeed eliminate low-mass 
halos through the suppression of short-wavelength fluctuations, as predicted by 
the linear perturbation theory. But to the contrary of expectation, our simula- 
tions yield singular cores in the collapsed halos, where the halo density profile i s 
similar, but not identical, to the NFW profile (jNavarro. Frenk fc Whitd 119971 ) . 
Such a profile arises regardless of whether the halo forms through accretion or 
merger. In addition, the virialized halos exhibit anisotropic turbulence inside a 
well-defined virial boundary. Much like the velocity dispersion of standard dark 
matter particles, turbulence is dominated by the random radial flow in most 
part of the halos and becomes isotropic toward the halo cores. Consequently the 
three-dimensional collapsed halo mass distribution can deviate from spherical 
symmetry, as the cold dark matter halo does. 
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Introduction 



Observations of low surface brightness ga laxies and dwarf galax i es indicate that the cores 
of ga lactic halos have shallow density profiles ( Dalcanton et al.lll997 



Swaters. M adore & Tre whella 



2000) instead of central cusps predicted by cold dark matter (CDM) (INavarro. Frenk fc White 
19971 ). In addition, the number density of dwarf galaxies in Local Gr oup turns out to b e 
an order of magnitude fewer than that produced by CDM simulations (IKlypin et al.lll999l ). 
These two features cast doubt on the validity of standard CDM. There have been at least 
three different solutions proposed to resolve these problems: (1) warm dark matter, (2) 
collisional dark matter and (3) fuzzy dark matter. 

Warm dark matter can suppress small-scale structures by free streaming. It seems to be 
able to both solve the over- abundance problem of dwarf galaxies and the singular core prob- 
lem. In this mode l the flat core is embedded within a radiu s a couple of percents of the virial 
radius (I Jingl l200ll ; IColins. Valenzuela &: A vila-Reesd 120081 ). and the core smoothly connects 
to the NFW profile (INavarro. Frenk &: Whitdll997l ) outside. Howe yer, this modification may 
gener ally adversely affect structures of somewhat larger scales (IHu. Barkana fc Gruzinov 
20001 ) . despite that fine tuning of the thermal velocity of dark matter p articles may still 



be able to have the larger scale structures consistent with observations (lAbazajianl 12006 
Viel et all 12008^ 



For collisional dark matter, the halo core can be flattened and dwarf gala xies destroyed, 
and N-body simulations confirm this conjecture (jSpergel fc Steinhardtll2000l ). But simula- 
tions also show that very freque nt collisions can yie ld even more singular cores than the 
standard collisionless CDM does (jYoshida et al.ll2000l ). This opposite behavior is indicative 
of the requirement of fine tuning for collisional parameters. 

The third solution to the problem is to treat dark matter as an extremely light bosonic 

dark matter (ELBDM) or fuzzy dark matter (IPress fc Rydenl|l990l ; ISinl|l994j ; lHu. Barkana fc Gruzinov 
20001 ). Axion has been thought to be a candidate of light bosonic dark matter. But for the 
light dark matter to erase the singular galactic core and suppresses low-mass halos, the 
particle mass must be far smaller than axion (m~ 10 -22 eV), so low that the uncertainty 
principle operates on the astronomical length scale. Much like axions, the ELBDM is in a 
Bose-Einstein condense state produced in the early universe. These extremely light particles 
share the common ground state and is described by a single coherent wave function. Its 
de-Broglie wavelength is comparable to or even somewhat smaller than the Jean's length 
( iDavies fc Widrowill997l ). where the quantum fluctuation provides effective pressure against 
self-gravity. 



Several previous works have pondered on such an idea or its variance (ISin 



1994: 



Hu. Barkana fc Gruzinc 
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2000l ; ISiddhartha Urena-Lopea 120031 ). in which the wave mechanics is described by the 
Schrodinger-Poisson equation with Newtonian gravity or by the Klein-Gordon equation with 
gravity. The Schrodinger-Poisson system addresses the scale-free regime of quantum me- 
chanics, where the Jean's length is a dynamical running parameter. On the other hand, 
the Klein-Gordon system makes use of the Compton wavelength as a natural length scale 
to create the flat core in a halo. Widrow & Kaiser conducted simulations for the two- 
dimensi onal Schrodinger-Poisson system to approximate the standard collisionless cold dark 
matter (IWidrow fc Kaiser! Il993l ). In the 2D case, the 1/r gravitational potential is replaced 
by log(r), and the 2D force law in their simulation becomes of longer range than it actually is 
in 3 D. Due to the lack of 3D numerical simulati ons, some au thors resort to spherical symme- 
try (jsin 1994 ; Siddhartha &: Urena-Lopez 2003 ) or even ID ( Hu. Barkana &: Gruzinov 2000 ) 
to study this problem. These simplifications may not capture what actually results in a 3D 
system with realistic initial conditions. In particular, the existence of a flattened core has 
been derived or inferred from these previous works of ID system or with spherical symmetry. 
In this paper we report high-resolution fully 3D simulations for this problem. Surprisingly, 
our simulations reveal that the singular cores of bound objects remain to exist even when 
the core size is much smaller than the Jean's length. 

In Sec. 2, we provide an explanation for the possible existence of the Bose-Einstein 
state for the extremely low mass bosons under investigation here. We then discuss two 
different representations of ELBDM and the evolution of linear perturbations for the two 
representations. In Sec. 3, the numerical scheme and initial condition are described. We 
present the simulation results in Sec. 4. In Sec. 5, we look into the physics of collapsed cores 
with detailed analyses from different perspectives. Finally the conclusion is given in Sec. 6. 
In the Appendix, we present results of ID and 2D simulations and demonstrate singular 
cores do not occur in ID and 2D cases. 



2. Theory 

2.1. Bose-Einstein Condensate 

A Bose-Einstein condensate (BEC) is a state of bosons cooled to a temperature below 
the critical temperature. BEC happens after a phase transition where a large fraction of 
particles condense into the ground state, at which point quantum effects, such as interference, 
become apparent on a macroscopic scale. The critical temperature for a ga s consisting of 
non- interacting relativistic particles is given by (IBurakovsky fc Horwitzlll996l ): 



K 3m J ' 



X 
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where the Boltzmann's constant and speed of light have been set to unity. Given the ex- 
tremely low particle mass assumed here, T c is derived from the relativitic Bose-Einstein 
particle-antiparticle distribution with the chemical potential set to particle mass m. Here 
the "charge" density n c h = n + — n_, where n + and n_ are the number densities of particles 
and antiparticles in excited states. On the other hand, we have n c h ~ (m/T)n + , and it 
follows that T c ~ (|^) 1//2 . Note that n + scales as a" 3 and T as a -1 , and it follows T c scales 
as a -1 . It means that when T is below T c at some time after a phase transition, the temper- 
ature will remain sub-critical in any later epoch. As an estimate, if we assume one percent 
of ELBDM to be in the excited states after its decoupling, the current critical temperature 
becomes 

T c = 3 x lO-^r^i^-y^eV. (2) 
ev ev 

Substituting m ~ 10~ 22 eV and T ~ 10~ 4 eV, the same as the present photon temperature, 
we find that the current critical temperature T c = 0.3eV » T. Hence ELBDM, if exists 
and accounts for the dark matter, may very well be in the BEC state ever since a phase 
transition in the early universe. Despite ELBDM particles in the excited state are with a 
relativistic temperature, almost all particles are in the ground state and described by a single 
non-relativistic wave function. 



2.2. Basic Analysis 

The Lagrangian of non-relativistic scalar field in the comoving frame is 

and t he equation of motion for this Lag rangian gives a modified form of Schrodinger's Equa- 
tion ( jSiddhartha fc Urena-L6pezll2003l ) : 

Oib ft? 

ih-^- = VV + mVi(), (4) 

at 2a 2 m 

where ip = 0(rio/a 3 ) _1/ ' 2 with </> being the ordinary wave function, no the present background 
number density and V is the self-gravitational potential obeying the Poisson equation, 

VV = 47cGa 2 5p = —p {\^\ 2 - 1). (5) 

a 

The only modification to the conventional Schrodinger-Poisson equation is the appearance 
cr 1 associated with the comoving spatial gradient V, and the probability density |?/>| 2 to be 
normalized to the background proper density p/m. In the above, 

3# 2 



Po 



8ttG 



= mn Q (6) 
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is the background mass density of the universe. 

To explore the nature of the ELBDM, we first adopt the hydrodynamical description 
to investigate its linear evolution. This approach is not only more intuitive than the wave 
function description, its advantage will also become apparent later. Let the wave function 
be 

V> = A /fe<f, (7) 
V no 

where n = na 3 , the comoving number density and n is the proper number density. The 
quadrature of Schrodinger's equation can be split into real and imaginary parts, which be- 
come the equations of accerleration and density separately, 

d 1 W ft 2 „ r V 2 ^, 

^ v+ ^ v " Vv + — -« v( ^r ) = (8) 

| + >>v)=0, (9) 

where v = VS/m is the fluid velocity. There is a new term depending on the third-order 
spatial derivative of the wave amplitude y/n in the otherwise cold-fluid force equation. This 
term results from the "quantum stress" that acts against gravity, an d it can be cast int o 



a stress tensor in the energy and momentum conservation equation (jChiuehl Il998l . 120001 ). 
The quantum stress becomes effective only when the spatial gradient of the structure is 
sufficiently large. 

The fluid equations, Eqs.(5),(8) and (9), are linearized and combined to yield 

e_ a ^ Sn _™^ Sn+ ^ m = (l (10) 

Upon spatially Fourier transforming 8n, it follows 

d 2 dn k 3H 2 n m h 2 k A 

-j; a —j2 ( — s + n k = 0, (11) 



dt dt 2a Am 2 a 



which can be recast into 



d 2 

x 2 —n k + (x 2 -6)n k = 0, (12) 

and the solution to this equation is 

(3 cos x — x 2 cos x + 3x sin x) 
n k = - 2 (13) 

where x = hk 2 / {m\J H 2 a) and a = (t/t ) 2 ^ and f2 m = 1 appropriate for early universe 
have been assumed. In the small-/c limit, x is small and n k ~ x~ 2 , which grows in time as 
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a; for large x the solution oscillates. Fig.(l) depicts the solution, Eq.(13). From Eq.(12) we 
can easily identify the oscillating solutions when x 2 > 6, thereby defining the Jeans wave 
number: 

kj = (6a)^(!^)V2. (M) 

Beyond the Jeans wave number, the perturbation is suppressed by quantum stress. 
Moreover, the Jeans wave number scales as a 1 / 4 and is proportional to m 1 / 2 (Hu, Barkana & 
Gruzinov 2000). We shall come back in a later section to examine up to what evolutionary 
stage the linear solution of n k can remain valid. 

Next, we linearize the Schrodinger-Poisson equation to derive a governing equation for 
an alternative wave-function representation. The wave function can be separated into real 
and imaginary parts, 

ip = l + R + U, (15) 

with fi, / <C 1. In the linear regime, we have the real part of linearized Schrodinger's 
equation 



and the imaginary part 

The Poisson Equation becomes 



S f = «V*I. (17) 

W = —p R. (18) 
a 



The spatial Fourier components of gravitational potential and ip satisfy 

8nGpoR k 
a k 



d h 2 k 2 

~ h tJ* = ^^ R * + mV *> ( 2 °) 
dt 2a z m 

and 

2a 2 mdR k 

h = -hk^^f (21) 
Combing the above, we have, as Eq. (11), that 

d 2 d „ ,3H 2 n msn , h 2 k 4 n . . 

X« 2 X^ - i—TTT 1 )^ + = 0, (22) 



and Rk has, up to a constant factor, the same solution as n k . Note that since dR k /dt = dR k /a 
for low-A; modes, it follows that \I k \ = (2mH a 1 / 2 /hk 2 )\R k \ = (k/ /y/3k?)\R k \ > \R k \. This 
feature will serve as one of the indicators for the validity of the linear regime in the wave 
function representation. 
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3. Numerical Scheme and Simulations 

3.1. Numerical Scheme 

We normalize the length to the computational grid size, A, and further define rj = 
(mA 2 H ) I h. The value of rj determines the size of Jean's length relative to the computational 
grid size. Set V = ^V. The dimensionless Schrodinger-Poisson equations becomes 

V = -2^ V ^ + ^T^' (23) 

and 

V 2 U = (H 2 -l), (24) 
where U(x) = V(x)/( 3 2 ™ v ), the dimensionless gravitational potential, and r = H t. 

Given a Hamiltonian H, one can evolve the wave function through (23). It is simply a 
unitary transformation of the system, 

= e~ mdt ij j . (25) 

We use the pseudo-spectral method to solve the Schrodinger equation in the comoving box. 
Let K be the kinetic energy operator (K = — ^V 2 — > k 2 /2rja 2 in Fourier space) and W 
the potential operator(W = in real space). The evolution is then split into 



e -awt = e -i(K+w)dt = i _ f(K + W)dt - -(K 2 + KW + WK + W 2 )dt 2 + 0(dt 3 ). (26) 

On the other hand, we need to consider the non-commutative relation between K and W, 
where 



e -iKdt e -iWdt = 1 _ t(K + w ^ dt _ l _¥?dt 2 - -W 2 dt 2 - KWdt 2 + 0{dt 3 ), (27) 

e -iWdt e ^Kdt = x _ ^ w + jQfc _ }_ W 2 df 2 _ i K 2 rft 2 _ WKdf 2 + Q(dt 3 ). (28) 

2 2 

It follows that we obtain, to the second order accuracy, 

e -i(K+W)dt ^ 1 ^ e -iKdt e -iWdt + e -iWdt e -iKdt^ ^9) 

2 

which will be adopted to advance the time steps. 
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For each time step, the kinetic energy operator is calculated in the Fourier domain, 

^■+1 = e ~^ dt ^ (30) 
and ip is advanced in real space with the potential energy operator, 

^(x) J+1 = e-^^W 7 - (31) 



To ensure numerical stability, we restrict the magnitude of dt that rotates the phase angle 
of wave function less than ~ in each time step, 

7r ra 2 , . 

dt < 777"^ — 2 ' ( 32 ) 



2 fe 



max 



7TO . . 

dt< — . 33 

~~ fifi nil 

I'm 1 1^ max 

In the early stage(a ~ 10~ 3 ), the stability condition is governed by the kinetic energy term, 
where k max 2 = 3tt 2 and dt < (Qn)^ 1 (rja 2 ) . At the late time, the gravitational potential 
becomes ever deeper, and therefore dt is controlled by the potential energy, where U max is 
the greatest value of potential in the real space. 



3.2. Simulation Scale 



We prepare the initial conditions with CMBFAST (ISeljak &: Za 



darriaga]ll9961) at z 



Hu. Barkana fc Gruzinov 



1000 w ith ACDM cosmology. Such initial conditions differ from that of 
(120001 ). where the Compton length of ELBDM already has imprints on the power spectrum 
at z = 1000. We choose this initial condition because only a few low- A; modes can grow for 
our choice of Jean's length and the details of initial power spectrum are irrelevant at the 
late time. The simulations run up to 1024 3 -grid resolution in a comoving box. 

For simulations in a much larger box, the background density averaged over an lh- 1 Mpc 
box can often change with time, the so-called environment effects, where galaxies prefer to 
form in regions of high background density. Here, we ignore the environment effect by fixing 
Q m = 0.3. We let the dimensionless parameter 77 = 1.22 x 10 -2 and 4.88 x 10 -2 for the 
1024 3 and 512 3 simulation boxes, which give a Jeans wavelength 50 kpc at z — 0. This value 
of 7] corresponds tom~ 2.5 x 10~ 22 eV. In the rest of this paper, we shall report only the 
simulation results of the highest resolution. 
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4. Results 

4.1. Validity of Linear Perturbation Theory 

We depict the early evolution of density power spectrum in Fig (2) from z = 100 to 
z = 10. The density power spectrum increases as a 2 for modes with wave number <C kj, 
and kj indeed increases as a 1//4 . These features are in agreement with what the linear 
perturbation theory predicts. 

On the other hand, we show the early evolution of 4|_R fc | 2 and 4|/ fe | 2 in Fig. (3) from 
z = 400 to z — 10. As expected the low- A; modes grow initially, while the high-A; modes 
are suppressed by quantum pressure. It is surprising to find that the linear evolution of 
ipk is valid only for a short period of time before z = 200. After that, the wave function 
deviates from what the linear theory predicts. In particular, the linear theory of wave 
function representation predicts that R k and I k grow as a and a 3//2 respectively, and \Ik\ 2 = 
(kj 4 /3k 4 )\Rk\ 2 ^> \Rk\ 2 for low-fc modes. This prediction agrees with the simulation result 
only when z > 200. At a somewhat later epoch than z = 200, we observe that the difference 
between \I k \ 2 and \R k \ 2 diminishes, and at z — 10 we find \I k \ 2 ~ |-Rfc| 2 - 

This problem is also manifested in the growth rate of the linear density power spectrum 
\n k \ 2 a: 4|i? fc | 2 . It is found that 4|i? fc | 2 indeed scales as a 2 before z = 200. When z < 200, the 
quantity 4|i4| 2 grows much faster than a 2 , and 4|i?fc| 2 becomes about one order of magnitude 
greater than what the linear theory predicts at z — 100. That is, the density power spectrum 
| ?7-fc | 2 vastly deviates from, and is much less than, 4|i? fc | 2 even since early on in the evolution. 

To examine this peculiar feature, we construct the imaginary part of wave function J(x) 
from I k of a few lowest k modes and depict |(/ 2 )fc| 2 on the same plot as 4|i? fc | 2 at z — 100 
in Fig. (3). It is found that |(/ 2 ) fc | 2 coincides with 4|_R fc | 2 at low-fc. Since n k ps 2R k + {I 2 ) k , 
it follows that (I 2 ) k has the opposite sign but approximately the same magnitude as 2R k so 
that the two terms of n k almost cancel to yield \n k \ 2 4|i?fc| 2 . Thus nonlinearity already 
sets in as early as z = 100 in the wave function representation. On the other hand, the 
fluid representation does not have such a problem. We find that \n k \ 2 of low- A; modes in the 
fluid representation agrees with what the linear fluid theory predicts even as late as z = 1, 
despite that the high-A; modes already become nonlinear. The difference arises from that 
(VS/m)(= v) in the fluid representation remains small at low- A;, in contrast to I — S in the 
wave function representation which has a large amplitude for low-fc modes even at high z. 
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4.2. Weakly Nonlinear ity Regime 

Shown in Fig. (4) is the evolution of |n fc | 2 (= \2R k + (R 2 + I 2 )k\ 2 ) for 1.5 < z < 10. The 
initial of high-A; mode has been linearly suppressed and is later replaced by the high- 
k modes that are nonlinearly generated beginning around z — 5. The nonlinear coupling 
arises from the coupling Vip in the Schrodinger equation. Since V is dominated by low- A; 
modes, the nonlinear coupling transfers modal energy local in k space from one mode to 
the neighboring mode, and from low k to high k. The gravitational potential V is barely 
evolving in the weakly nonlinear regime, and hence the dynamics in this regime is for the 
wave function to settle into almost static potential wells. 

We note that if V were exactly static, the Schrodinger equation would have been linear 
and the coupling from low to high k modes would simply be the linear evolution of wave 
function starting with a non-eigenstate. This argument may explain why the linear theory 
of low- A; modes works so well even after high-A; modes are excited in this regime. Shown 
in the left panel of Fig. (5) is the real-space configuration of Sn at z = 3 in this weakly 
nonlinear regime, where the wave function settling into individual quasi-static potentials 
well is underway. 

4.3. Strong Nonlinearity Regime 

Plotted in Fig. (6) is the evolution of density power spectrum after z = 1.5. This is 
the strongly nonlinear regime where the gravitational potential develops deep wells at the 
collapsed cores. To illustrate the contribution of the few collapsed objects to the final high-A; 
power spectrum, P(k,0), at z — 0, we remove all matter outside the virial radii of these 
collapsed halos and construct the power spectrum of these artificial objects, P h (k,0). The 
power spectrum of the removed matter, P b (k,0), is also constructed for reference. The two 
power spectra along with the original power spectrum are depicted in Fig. (7) for comparison. 
Clearly the bound objects contribute to almost all the power contained in the original power 
spectrum, except for the low- A; modes that are contributed dominantly by Pb(k, 0). These 
few lowest-A; modes are grown out of the initial noise and remain so in the final configuration. 
That is, despite that the initial condition possesses many independent degrees of freedom, 
the final configuration has only a few degrees of freedom, where almost all randomly placed, 
small collapsed objects seen in standard CDM simulations are entirely suppressed. 

During the final collapse phase, the core undergoes large-scale oscillations that send 
out waves to remove excess angular momentum deposited into the core region, rendering 
the core to settle into an almost stationary configuration in the physical space. Fig. (8) 
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shows the wavy structures of this nature around the collapsed halo A at z=l. Even in the 
quasi- stationary state of these halos at z — 0, we find that this wave phenomenon is still 
pronounced around the halos, as will be discussed later. 

There are two collapsed halos, A and B, of mass 5.7x 10 9 Mq and 5x 10 9 Mq respectively, 
at z = in our simulation as shown in the right panel of Fig. (5). Halo B is subject to major 
merger around z = 0.7. Shown in Fig. (9) are halo B before and after the major merger. The 
final density profiles of halos A and B are plotted in Fig. (10). Interestingly, they both develop 
singular cores, in spite of the presence of quantum pressure. Both power-law singular cores 



have a power index —1.4, reminiscent of that of the standard cold dark matter (IGhigna et al. 



20001 ). The density profiles of the outskirts also obey power law, with a power index —2.5, 



slight ly shallower than that produced by the cold dark matter (jNavarro. Frenk fc White 



19971 ). Similar density profiles arising from different formation processes, i.e., accretion 



versus merger, suggest that the profile can be universal. 

Note that the final collapsed core contains angular momentum through angular depen- 
dence in the wave function. The angular dependence manifests itself with large-amplitude, 
small-scale fluctuations in the wave function. To examine this aspect of the halo, we let the 
wave function be represented by ip = fe tS . The specific kinetic energy is obtained through 
the real part of the expression: 

^ u (ysy v 2 / . v/ VQ v 2 s 

- = ~^r ] ~ ( t ' vs + ~2~ )] - (34) 

On the other hand, the specific flow energy can be evaluated through the real part of the 
following: 

V(*V|*|>) _ 1 r (VS)' ,35) 



47? 2 (V>7H 2 )) V 2 2 v 2 
Combining the two, the specific internal energy is obtained. Plotted in Figs. (11) and 
(13) are the 2D slices of density for halos A and B, showing large-amplitude, small-scale 
fluctuations in density. We also plot the same slice for the 2D flow velocity V±S/r] (= 
— iV ±{il) 2 /\^\ 2 ) / {2r]ip 2 /\ip\ 2 ) in Figs. (12) and (14). The velocity patterns clearly reveal well- 
defined boundaries of turbulence regions in halos A and B against the infall. The flow 
becomes randomized inside this sharp boundary. Despite the boundary outlines an accre- 
tion shock-like structure, we find in the density slice that there is no obvious jump at the 
boundary and it is not a shock. Thus, there is no analogy of such a structure in the fluid 
system. This peculiar feature warrants our further investigation in the future. 

We next investigate the virial conditions in the two halos. Plotted in Figs. (15a) and 
(15b) are the ratios of the kinetic energy integrated up to a radius r to the potential energy 
(Jq Anr' 2 {?>VL m rj/4)[n{U — U min )}dr') integrated up to r for halos A and B. The virial ratios 
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are 0.5 at the average radii of the infall boundaries, within which turbulence occurs. In 
addition, we also plot the specific flow energy and the specific internal energy respectively 
in Fig. (15). It is found that the specific internal energy is twice as large as the specific flow 
kinetic energy in the interiors of two halos. 

Virialization can be correlated with the flow equi-partition. Plotted in Fig. (16) are the 
random tangential specific flow energy and the random radial specific flow (subtracted off 
the mean radial infall) energy averaged over spherical shells for the two halos. The tangential 
flow energy is about twice the radial flow energy only well within the halos, thus providing 
an evidence of equi-partition at the halo cores. In the outskirts of the halos, the random 
radial flow energy is larger than the equi-partition value. This aspect is reminiscent of the 
velocity dispersion in a standard CDM halo. 

Equi-partition is also related to the sphericity of mass distribution. Significant large- 
scale angular dependence in the wave function can yield aspherical halos. We define the 
quadrapole-to-monopole ratio as 



n = (A 1 -A 2 ) 2 + (A 2 -A 3 ) 2 + (A 3 -A 1 ) 2 
2(A 2 + A 2 + A 2 ) } 



(36) 



where Ai, A 2 and A3 are the eigenvalues of J (rr/r ,3 )nd 3 r, with (3 = 3.5 to weigh in favor 
of the core. The Q value characterizes the low-order angular dependence of wave func- 
tion, which assumes the extreme value zero or unity when the density profile has spher- 
ical symmetry or a one- dimensional shape. It is found that Q = 0.34 for halo A and 
0.11 for halo B respectively. Perhaps violent relaxation after a major merger acceler- 
ates halo B to assume spherical symmetry. We note that the quan tum stress is in fact 
anisotropi c: T?- = [d^/n) {d js /n)/rf - <%(V 2 n)/4r/ 2 Jchiuehl boooh . Unlike fluid dark 



matter (jYoshida et al.l |2000| ) , the density asphericity arises from the anisotropic stress of 
quantum mechanics, similar to that produced by collisionless dark matter. The similarity 
between the quan t um d ynamics and the collisionless particle dynamics in fact motivated 
Widrow Kaiserl (119931 ) to propose a model that approximates the latter by the former. 



5. Conclusion 

As far as we know of, this work presents the first result for the stud y of Bose- Einstein con- 
densa te under self-gravity via high resolution (1024 3 grids) simulation. iHu. Barkana fc Gruzinov 



( 120001 ) conjectured that if the dark matter is ELBDM, it can solve the long standing problem 
of far too many low-mass halos present in the standard CDM simulations, and also explains 
the existence of flat cores in some galaxies. In this work, we confirm that low mass halos are 
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indeed suppressed by quantum stress even when the small-scale fluctuations are abundant 
in the initial power spectrum. This result is a consequence of long-time linear suppression of 
the small-scale modes. We also find, from our simulations of different grid resolutions, that 
collapsed halos develop singular cores regardless of the halo formation processes. All these 
runs produce convergent density profiles. Our 1024 3 highest resolution run gives singular 
density profiles similar to what standard CDM simulations produce. 

In retrospect this singular-core result may not be too surprising, as it arises from an 
almost scale-free Schrodinger-Poisson system. This system is not exactly scale free because 
there exists a Jeans length for fluctuations that are small in amplitude. However, when the 
local density much exceeds the background density, the latter becomes locally ill-defined, 
the Jeans length no longer has any physical significance, and the Schrodinger-Poisson system 
becomes locally scale free. Being locally scale-free, the system develops singularities within a 
finite time. By contrast, the conservation of phase-space density in classical particle dynamics 
precludes the space density of standard dark matt er particles from developing any singularity 
( jChiueh fc Wool 119971 ; iDalcanton fc Hogan Il200ll ). and explains the existence of a flat core 
in the warm dark matter model. Note that such a phase-space constraint does not exist for 
nonlinear wave dynamics; one example of this nature is a system described by the non linear 
Schrodinger equation with attractive self- interaction (jSulem. C. fc Sulem. P.-L.lll999l ). 



Most recent observations of rotation curve in low-surface-brightness galaxies indicate in- 
conclusive results, as far as the existence of singular halo core is concerned. Some galaxies are 
claimed not to possess singular halo cores, and some are if non-circular motion is taken into 



with the constraint given by ACDM cosmology ( 


Swaters et al. 


200 


3: 


Zakurisson et al. 


2006; 


Kuzio de Narav et al. 


2006; 


Kuzio de Narav, McGaugh & de Blok 


2008 


). Given the present 



status of observations, if galaxies indeed contain singular cores, ELBDM will likely be the 
only viable candidate for the dark matter that, on one hand, permits the galactic-scale, 
NFW-like halo cores, and on the other hand suppresses the sub-galactic low-mass halos. 
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APPENDIX 



A. Comparison with One and Two Dimensional Simulations: 

We also examine the one- dimensional and two dimensional simulations with the same 
numerical scheme used for the three-dimensional simulations. It is hoped that we can verify 
the scheme by reproducing the results of previous works. 

So far as we know of, the only work on one-d i mensi onal simulation of the same physical 



system is the work of iHu. Barkana fc Gruzinovl (120001 ) . In this work, they simulated the 



ID sheet collapse in a three-dimensional expanding universe. Two cases of different Jeans 
lengths were examined. For the short Jeans length case, they found fast oscillating solutions, 
in time and space, which corresponds to a solution of repeated collapses and rebounds. For 
the long Jeans length case, they found slow oscillating solution. In either case, the collapsed 
sheet cores are smooth with no cusp. 



The only two-dimensional simulation in the literature is that of lWidrow fc Kaiser! (119931 ). 
Unfortunately, they did not examine individual profiles of collapsed filaments and no com- 
parison can be made, since the focus of this work placed on the equivalence of the light 
bosonic system and the cold dark matter system. 



This preprint was prepared with the AAS IATj^X macros v5.2. 
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We run one-dimensional simulation with 8192 grids and a small-amplitude initial random 
fluctuation. The evolution follows the linear growth described in the text. Nonlinearity 
sets in quite quickly when the density fluctuation grows to a amplitude comparable to the 
background density. The saturation amplitude of Sp/p is around 10. The solution indeed 
oscillates, corresponding to repeated collapses and rebounds. We show the solutions in 
Fig.(Al) for two nearby time steps where Aa = 0.001 near a = 1. It is apparent that the 
solution has no singular c ore. This solution is similar to what was found in the work of 



Hu. Barkana fc Gruzinovl (120001 ) . 



The two-dimensional case has also been examined with the same scheme as the 3D code 
with 4096 2 grids. We find that the collapsed filaments at z = do not contain singular cores; 
rather, several dense clumps of Sp/p ~ 1000 are scattered over a area of finite radius in the 
core region, as shown in the zoom-in image in Fig.(A2). 

It is understandable why low dimensional objects do not develop singularities. This 
is due primarily to weakening of the focusing power in a three-dimensionally expanding 
universe. It appears that 2D is the critical dimension, where the singularities are just not to 
appear. In 3D, the self-focusing power can be so strong as to bind the dense clumps tightly 
to form singularity. 
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Fig. 1. — The solution of the linear perturbation given by Eq.(6). The vertical line labels 
the location of squared scaled Jean's wavenumber at x — 6. 
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Fig. 2. — The linear evolution of the power spectrum from z = 100 to z = 10 in 
box. The low-fc power obeys the linear scaling oc a 2 . 



an 1 A 1 Mpc 



-19- 



4R 2 (400) 



9--o 4f (400) 
— — 4R 2 (200) 
4? (200) 




s\ , 1 # 1 1 

J( U001 0.01 0.1 1 10 

k (h kpc) 

Fig. 3. — Evolution of 4\Rk\ 2 (z) and 4|/fc| 2 (,2) to check against the linear theory. Deviation 
from the linear theory is evident from z=200 on. Black dots are |(/ 2 )fc| 2 at z = 100 con- 
structed from few low-fc modes to show the cancellation between 2R^ and (I 2 )k so as to make 
n k < 2R k . 
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Fig. 4. — The weakly nonlinear evolution of the power spectrum from z = 10 to z = 1.5, 
where the high-fc modes are seen to be nonlinearly excited. 





Fig. 5. — Two-dimensional projections of density in real space in a lh 1 Mpc comoving box 
at z=3 (left panel) and z=0 (right panel). Halo A and B are at the top left and bottom left. 
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Fig. 6. — The nonlinear evolution of the power spectrum from z = 1.5 to z — 0.0. The 
highest-/c modes acquire their full power after z = 0.4, indicative of the creation of singular 
halo cores. 
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Fig. 7. — The comparison of the halo power spectrum Ph (circle), the background power 
spectrum P b (triangle) and the full power spectrum P (square) at z = 0. Note that Ph 
matches P at high k and P& matches P at low k. 




Fig. 8. — Waves are sent out from the collapsing halo A at z — 1 (left panel); it develops an 
oblate singular halo at z = (right panel). 
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Fig. 9. — Halo B is subject to major merger at z = 0.7. The left panel reveals the progenitors 
at z = 1 and the right panel shows a singular halo with high degree of spherical symmetry 
at z = 0. 
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Fig. 10. — Density profiles of two massive halos at z=0. The left panel plots the profile of 
halo A and the right panel the profile of halo B. In both panels dot-dash lines and solid lines 
denote the power law of indices —1.4 and —2.5, respectively. 
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Fig. 11. — A two-dimensional slice of density for halo A through the core. 
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Fig. 13. — A two-dimensional slice of density for halo B through the core 
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Fig. 15. — Virial ratios of the kinetic energy integrated up to a radius r (square) to the 
potential energy integrated up to r, the specific flow energy (plus) and specific internal 
energy (star) for halos A (upper panel) and B (lower panel). This virial ratios are about 0.5 
at the average radii of the infall boundaries. The specific internal energy is about twice as 
large as the specific flow kinetic energy in the interiors of the two halos. 
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Fig. 16. — Specific radial flow energy and tangential flow energy of halo A (upper panel) 
and halo B (lower panel). 
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Fig. A.l. — The density profiles of two nearby time steps where Aa = 0.001 near a = 1 in a 
ID simulation of 8192 grids. 
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